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ABSTRACT 

In this paper we re-analyse the amplification process of broadband continuum radiation by 
astronomical masers in one-dimensional case. The basic equations appropriate for the scalar 
maser and the random nature of the maser radiation field are derived from basic physical 
principles. Comparision with the standard radiation transfer equation allows us to examine 
the underlying assumptions involved in the current theory of astronomical masers. Simula- 
tions are carried out to follow the amplification of different realisations of the broadband 
background radiation by the maser. The observable quantities such as intensity, spectral line 
profile are obtained by averaging over an ensemble of the emerging radiation correspond- 
ing to the amplified background radiation field. Our simulations show that the fluctuations of 
the radiation field inside the astronomical maser deviates significantly from Gaussian statis- 
tics even when the maser is only partially saturated. Coupling between different frequency 
modes and the population pulsing are shown to have increasing importance in the transport 
of maser radiation as the maser approaches saturation. Our results suggest that the standard 
formulation of radiation transfer provides a satisfactory description of the intensity and the 
line narrowing effect in the unsaturated and partially saturated masers within the framework 
of one-dimensional model. Howerver, the application of the same formulation to the strong 
saturation regime should be considered with caution. 
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1 INTRODUCTION 

Since its discovery in the interstellar medium, maser emission has been used as probe to study the different aspects of astronomical en- 
vironments. The strong maser emission has even provided detailed pictures of the central region of distant galaxies. There had been great 
expectation that masers could also provide information on the physical properties of the molecular gas at scales inaccessible by other obser- 
vational means. Theoretical studies have been carried out to find the pumping mecanisms and the necessary conditions for the existence of 
masers. However, due to the strongly nonlinear nature of masers, the relation between observed characteristics and the physical conditions 
of the masing medium is difficult to determine. 

At the fundamental level, the physics of astronomical maser is very similar to that of laboratory lasers, which are masers operating at the 
optical wavelengths. Most of the properties of laboratory lasers are well described by the semi-classical model of Lamb (1964). The feature 
that distinguishes laser radiation from classical light sources is the remarkable degree of spatial and temporal coherence largely due to the 
nature of stimulated emission and the use of high-Q cavity and other means of stabilization. The statistics of laser radiation has been shown 
to follow Poisson-like distribution (Arechi 1965). The similarity between the laboratory laser and the astronomical maser has prompted spec- 
ulations that maser radiation might not be Gaussian. However, previous observational attempts (Evans et al. 1972, Moran 1981) to measure 
the non-Gaussian statistics of maser radiation have proved negative. 

It was recognized very early in the study of astronomical masers (Litvak 1970) that the maser radiation field is broadband. The coherence 
time set by the bandwidth of the maser radiation is therefore much shorter than the typical time scale of molecular response. In deriving 
the transfer equations for maser radiation, two important assumptions are used by Litvak (1970) and later by Goldreich et al. (1973) in their 
study of the maser polarization. The first assumption specifies that the field is stationary and different spectral components are not correlated. 
The second assumption concerns the Gaussian distribution of the electric field components. Goldreich et al. (1973) reason that due to the 
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large dimension of astronomical masers, the phase shift between waves propagating in slightly different directions is so enormous that the 
waves become completely independent and random. The electric field seen by masing molecules is the superposition of independent and ran- 
dom electromagnetic waves and therefore will have Gaussian statistics according to the well-known central limit theorem. These two main 
assumptions allow considerable simplification in deriving the transfer equations. Although several theoretical studies of maser radiation are 
based on these assumptions, their validity has not been fully examined, especially in the saturation regime where analytical and numerical 
studies have been regularly carried out. We note here that the work by Field & Richardson (1984) and Field & Gray (1988) did try to take 
into account the partial coherence of the maser radiation in an approximate way. Their analysis suggests some difference in the predicted 
maser intensity in comparison with the commonly used maser theory. Furthermore, the recent debate on the theoretical descriptions of the 
polarization properties of astronomical masers (Watson 1994, Elitzur 1995) also underlines the urgent need to re-examine the foundations of 
the current maser theory. 

To simulate from first principles the amplification of radiation by a masing medium is quite complicated. The radiation field has to be con- 
sidered as an ensemble of all possible statistically independent realizations. Studying the amplification of maser radiation then requires a 
dramatic increase in computing power to follow the evolution of a large number of realizations of the incident radiation. However, from 
these samples we can study the effect of amplication on the statistical distribution of the radiation field and check the validity of the above- 
mentioned assumptions. The physical model to describe the amplification of broadband radiation by a one-dimensional maser has been 
formulated by Menegozzi & Lamb (1978). Using simple analysis they show that the Gaussian statistics can not be maintained as the radia- 
tion is amplified and begins to affect the population inversion of the maser. Menegozzi & Lamb carried out the simulations in several cases 
including the medium with Dopller broadenning, which is relevent to astronomical masers. Due to limited computing power, they could 
simulate only a handful of realizations of the radiation field, and therefore, the results are still inconclusive. It is surprising that no further 
study has followed the work of Menegozzi & Lamb (1978), particularly in view of the great computing power available today. 
A recent work by Gray & Bewley (2003) also deals with the broadband nature of the maser radiation field. The conclusion reached in their 
work indicates that different frequency components across the maser line profile are phase-locked even in the partial saturation regime. The 
phase-locking (or mode-locking) effect also means that maser radiation is pulsating. The result is unexpected because the phase-locking of 
different frequency modes was not seen in the work of Menegozzi & Lamb (1978). Furthermore, in laboratory laser systems the mode-locking 
can only be achieved through external control either with opto-electronic components (Q-switch) or with passive elements such as saturable 
absorbers. We could not easily find a clear explanation on the physical mechanism giving rise to this behaviour. We will discuss the model 
of Gray & Bewley and compare their result with that obtained from our simulations in the following sections. 

In this paper we will use the formulation worked out by Menegozzi & Lamb (1978) to study in detail the astronomical maser and simulate 
from first principles the amplification of broadband background radiation field by the maser. The characteristics of the maser radiation will 
be determined directly without recourse to any approximations. The generalization of our model to include the vector nature of the radiation 
field to study the polarization properties of masers is straightforward and will be described in a subsequent paper. We hope that our work will 
provide a better understanding of the underlying physical principles and the limitations of the current theory of astronomical masers. 

2 BASIC THEORY 

2.1 Statistical properties of radiation field 

The continuum radiation which propagates along the z axis and arrives on the boundary of the masing medium is the superposition of waves 
coming from many small, independent sources and therefore, is a random function of time. We assume that the incident radiation field E(t) is 
stationary, which means that the statistical properties of the field does not change with time. We further assume that the underlying probability 
distribution of random variable E{t) at any instant of time is described by a Gaussian with zero mean and variance < \E(t)\ 2 >=a 2 , where 
the brackets denote the ensemble average, or the average taken over a large number of realizations of the field. In practice we can measure 
the power of the field by averaging the output of a square-law detector: 



J -T/2 

Because E(t) is a random process, It defined as above for each realization of the field during the interval [—T/2, T/2] is also a random 
variable. As T tends to infinity, It will approach the average power of the radiation field: 



To relate the statistical properties to time average quantities we need to restrict the random process E(t) to be ergodic, i.e ensemble average 
is equal to time average of any member of the ensemble. Justifications for this assumption constitute the famous ergodic theorem (Doob 
1953). The ergodic property is also the underlying assumption of almost all observations in radio astronomy. Therefore we can write: 
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where the brackets denote the ensemble average. In the following sections of the paper we shall use the concept of spectral components of 
the random process E(t). For that purpose we might formally use Fourier transform to decompose E(t) into harmonic elements E(uj): 

/ + oo 
E(u))e- iat du 
-oo 

1 f + °° 

E(u) = ^ E{t)e-^dt (4) 

J — oo 

It becomes apparent immediately that because the random function E(t) does not decay to zero at infinity and the above integrals do not 
exist. In other words, the function E(t) is not absolutely integrable. We can circumvent this difficulty by truncating the function E(t) on a 
finite interval [—T/2, T/2] and define a new function Er(t) which is equal to E(t) over -T/2^ t and periodic outside this interval. 

We can now expand Er(t) into harmonic components using Fourier series: 



n— + oo 



E T (t) = J2 B K) e_i ""' (5) 

71.— — OO 

Where the angular frequency & 7l =2im/T and 

E{u n ) = i f ' E{t)e^dt (6) 



Because E(t) is a random function, the coefficients E(ui„) are new complex random variables. In general, the coefficients E(uj n ) are 
independent of each other only in the limit of large T, i.e < E(u>„)E(u} m )* >= for n 7^ m (Root & Pitcher 1955). For random process 
with continuum spectrum or "white noise", the values of E(t) at different instants of time are completely uncorrelated, or the auto-correlation 
is a delta function. Therefore the coefficients E(uj n ) are Gaussian random variables and independent even for finite values of the interval T 
(Goldman 1953, Rice 1944, Root & Pitcher 1955). Using Perceval's theorem we obtain the following relation: 

" =+° ° , pT/2 

\E{u> n )\ 2 = - / \E(t)\ 2 dt (7) 



T/2 



Physically, the term \E(uj n )\ 2 can be interpreted as the power of E(t) contained within the frequency interval Aw=w n +i —Ld n =2ir/T during 
the interval [—T/2, T/2]. By taking the ensemble average, which means averaging over a large number of realizations of the radiation field, 
we can define: 

^- < \E(lu„)\ 2 >= 5tK)Aw (8) 
Sir 

where c is the speed of light. In the limit T — » 00, the quantity St (w) approaches the power spectrum of the random process E(t) : 
lira SVH = S(w) (9) 

T — >oo 

The summation in the Fourier series can be formally replaced by the Fourier-Stieltjes integral in the limit that T tends to infinity by substi- 
tuting a n by dZ(uj), the increment of an orthogonal random process Z(uj), 

/+00 
e~ i[Jt dZ{ui) (10) 
- OO 

where the limit should be understood in the mean square sense, namely: 

/+00 
e~ iu>t dZ{u))f >= (11) 
- 00 

If the random function Z(u) is differentiable, i.e dZ(uS)/dw exists, the expression above will reduce to the normal Fourier transform. As we 
can see from Eq.[8] the magnitude of jdZ(tj)| is 0(V du>), many orders of magnitude larger than duj. Hence, the derivative of Z(ui) does 
not exist. More rigorous derivation of the spectral resprentation theorem for random processes can be found in standard texts such as Doob 
(1953), Yaglom (1962) and Priestley (1981). 

In the case of an emission line, the radiation is bandwidth limited around the frequency of the transition. If we take oj a b=27r^ a b, for the sake 
of simplicity, to be the central frequency of the maser line, the electric field of the radiation propagating along the 2 axis can be expressed as 
follows: 

E(z, t) = i {E(z, t)e- iu '«>«-" e > + c.c} (12) 

where c is the speed of light and the complex amplitudes E(z, t) are slow-varying functions of both time and space. The complex function 
E(z, t)e~' UJab ^ t ~ z ' c ^ is called the analytic signal representation of the electric field since the function is analytic in the upper half of the 
complex lu plane. In practice we can not simulate or measure the random process E(z,t) over an infinitely long interval of time, we shall 
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restrict ourselves to study realizations of the radiation field during an interval T long enough that provides adequate spectral resolution 
Alu=2tt/T. Using the spectral representation theorem, we can expand the time dependent electric field in terms of Fourier series: 

+00 

E{z,t)e- l ^ (t - z/c) = e — vtt'-z/c) ^ E(z,uj n )e-^ n{t - z/c) (13) 

n= — 00 

where the Fourier components are given by 

E(z,u> n )e iUn * /c = (1/T) f 7 {£(2,t)e- l ^ (t - z/c) }e^" t dt (14) 

J-T/2 

with u n = 27rn/T. For bandwidth limited signal the Fourier components E(z, u n ) will be zero for sufficiently large values of n. 

The complex components E(z,ui n ) are random variables, whose real and imaginary parts (E r , E 1 ) are jointly independent Gaussian with 

zero mean and the same variance < \E(z, uj n )\ 2 >/2. The joint probability distribution is (Mandel & Wolf 1965): 

P[E(z,u n )] = (1/tt < \E{z,lu„)\ 2 >)e X p[-(E r (z,iu n ) 2 + E\z,Lu n ) 2 )/ < \E{z,u n )\ 2 >] (15) 

Noting that the expression for the area element d 2 E = d(Re E)d(\m E) = ^d\E\ 2 d(f> in the complex plane of E(z,u n ) we can see that 
the joint density distribution for amplitude and phase of E(z,u n ) as given by Menegozzi & Lamb (1978) follows directly from the above 
expression: 

$[E(z,u n )] = (1/2tt < \E{z,uj n )\ 2 >)exp[-\E(z,u n )\ 2 / < \E(z,lu„)\ 2 >] (16) 
with the following normalization: 

<b[E{z,Lo n )]d\E{z,uj n )\ 2 = 1 (17) 

The complex components E(io n ) for each realization of the input continuum radiation field at the boundary of the masing medium can be 
generated with a random number generator. As described in Menegozzi & Lamb (1978), if x is a real random number uniformly distributed 
in the interval between and 1 (excluding the end points), the phase (f> is simply 2nx. The amplitude of the harmonic component is obtained 
by the relation: 

\E{z,uj n )\ 2 = - < \E{z,uj n )\ 2 > ln(l - x) (18) 

where < \E(z, ui n )\ 2 > is the variance or the intensity of the radiation field at frequency uj n and must be specified beforehand. 
In the studies of Goldreich et al. (1973) and Deguchi & Watson (1990) the electric field is assumed to be stationary and can be decomposed 
into harmonic components using the usual Fourier tranform in similar fashion to Eq. [4] The intensity or Stokes parameters, if we consider 
the vector nature of the radiation field, are defined as the ensemble averaged quantities involving harmonic components at the corresponding 
frequency, similar to our definition above. Although Eq.[4]has only symbolic meaning, it can be reinterpreted using truncation technique as 
presented above without any change in the subsequent calculations. 

The procedure used by Elitzur (1992) is markedly different because only the decomposition of a single realisation of the radition field into 
Fourier components is considered. By truncating the function E(t) on a finite interval [— T, T] and performing the usual Fourier transform, 
the function Ft{w) is defined as: 

F T (u) = J E(t)e iuJt dt (19) 
The power spectrum S(u>) of the field is defined as (Elitzur 1992, Eq. 4.1.7, 9 and 10): 

S(u)du= lim — / E{t) 2 dt= lim - / \F T (ui)\ 2 du (20) 

J — 1 J — 00 

The last equality follows from the Perceval's theorem for the truncated function E(t) over any finite interval T. The power spectrum S(u>) of 
the radiation field is then identified with the limit: 

S(u) = Km^ ^\Ft(lu)\ 2 (21) 

However, the existence of the average power of the stationary radiation field, i.e the limit of ^ J E(t) 2 dt as T tends to infinity, does not 
guarantee that the limit in the above equation exists. In fact, the value of 27r \Ft{uj)\ 2 /T, which is called the periodogram by its originator 
Schuster, fluctuates wildly around the true power spectrum S(cu). The periodogram is usually refered as the non-consistent estimator of 
the power spectrum in the signal processing literature. One way to achieve convergence is to perform a local averaging of the periodogram 
(Champeney 1989): 

l im f l im 1 f ^^IFtHI 2 ) = S(u) (22) 



On the theory of astronomical maser. I. Statistics ofmaser radiation 5 



This theorem shows that in order to obtain good estimate of the power spectrum, a trade-off in spectral resolution must be made, which is of 
course the principle of many smoothing schemes to estimate the power spectrum of random processes. 

If we want to properly define the power spectrum using time average of a single fluctuating function E(t), we need to resort to the techniques 
of generalized harmonic analysis (Wiener 1930). Simply speaking, to obtain the power spectrum we need to measure the auto-correlation 
function: 

1 f T/2 

R(t) = lira - / E(t)E*(t + r)dt (23) 

T^oo T J_ T/2 

The Wiener- Khintchin theorem then guarantees that the function S(lo), which can be identified with the power spectrum, exists and is related 
to the auto-correlation function R(t) by the usual Fourier transform: 

r+oo 

S(u)= / R(r)e-^ T dT (24) 



This procedure serves as basic principle of the XF or lag correlator (Thompson et al. 2001), which is currently more practical than the FX 
design and has widespread use in radio interferometry. 



2.2 Matter-radiation interaction in masers 

Since our aim is to simulate the spectral properties of the maser radiation, it is preferable to work from the beginning in the frequency 
domain. The basic equations are derived here for the idealized case of a scalar one-dimensional maser. That means the masing molecules are 
idealized as a system with only two (upper and lower) energy levels and the properties of the maser depends only on the z coordinate, which 
is also the propagation axis of the maser radiation. These assumptions are commonly used in prevous works on the theory of astronomical 
masers. Following Menegozzi & Lamb (1978), Goldreich et al. (1973), Deguchi & Watson (1990) we will derive the transfer equation for 
the radiation field within the framework of the rotating-wave approximation. That means when dealing with the frequency dependence of 
the electric field and polarization vector, we retain only the positive frequencies. The anti-resonant (negative) frequencies are ignored. As in 
previous section, we express the electric field and the induced polarization vector as: 

E(z,t) = t)e- < "-»<«-* /c >+c.c} (25) 

and 

POM) = i{P(z,t)e-^ ( '- z/c) + c.c} (26) 

Where uj ab is the center of the maser line, corresponding to the rest frequency of the transition between the upper level a and the ground 
level b of the line. The amplitude of the field E(z, t) and polarization P(z, i) are assumed to vary slowly with time in comparison to the term 
exp(-iujabt). As shown in the previous section, the amplitude of electric field and polarization vector can be expressed in terms of Fourier 
series for a time interval T, which corresponds to a spectral resolution of Auj = 2n/T: 

n — + oo 

E(z,t) = £(^n)e-^" (t - z/c) (27) 

n— — oo 

and 

n — + oo 

P(z,t) = P{z,uj n )e-^ {t - z/c) (28) 

n— — oo 

The radiation transfer equation (Goldreich et al. 1973, Deguchi & Watson 1990) can be written as follows: 

GM) *<*•<> =^.'> 

or in the frequency domain: 

— E(z,u„) = P{z,Lo n ) (30) 

dz c 

The polarization P describes how the masing medium interacts with the radiation field at the frequency of the maser line. To calculate its 
value we need to use the density matrix, which describes the masing medium. The density matrix at any position z for a group of molecules 
moving with velocity v can be written as: 

P(z,v,t)=( Paa Pab ) (31) 

V Pba Pbb J 

The Hermitian property of the density matrix implies that pba = Pab- When there is no possible confusion we will drop the explicit 
dependence of the density matrix and the radiation field on the spatial coordinate z. To work in the frequency domain, we need to expand the 
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density matrix elements into Fourier series, in a similar fashion to the electric field: 

n— + 00 

Pab(z,v,t) = e~ Mab ( *-*/°) [p ab (2,w,a; n )e^'* 2/c ]e"' u '' ,t 



p^(z,V,t) = \J Paa(2, u, w„)e 8Wn(< z/c) (32) 

n=— 00 

71 — + OO 

p bb (z,v,t) = 2J Pbb(z,w,tJ„)e" ! ^' l(t ~ z/c) 

n=— 00 

Because the diagonal elements of the density matrix p aa and pbb are real, we have the following relations p aa (2, v, ui)=p a& {z, v, — lj)* . The 
Hamiltonian of the molecules is the sum of the Hamiltonian Ho, representing isolated molecules, and the matter-radiation field interaction 
matrix V. The well known evolution equation of the density matrix p(z,v,t) (Goldreich et al. 1973, Deguchi & Watson 1990) has been 
shown to have the following compact form by Icsevgi & Lamb (1969) and Sargent et al. (1974): 



fd d 
.di +V d- 2 

d d \ , . , „ i 



{.W: + V lfz) = ~( iUJab + HPab - ^-V ab (pbb - Pa 



\8t + V ~dz ) = ~ ^ ( ° aa _ ^(Vab/Oba - V ba p a b) (33) 

— Jpti, = Ab(l> ) — Tb pbb — ^-(VbaPab — V a bPba) 

Where A a and Ab are the pumping rates into the upper and lower levels, respectively. r a and Yb are the loss rate due to pumping and 
collisional decoherence (Sargent et al. 1974). V a b are the matrix elements of the interaction matrix V. For the sake of simplicity, we assume 
here that the loss rates are the same for the lower and upper levels. Substituting the Fourier expansion of density matrix elements into the 
above equations (except for the terms involving the interaction matrix V to be written explicitly later) and collecting term by term, we obtain: 

Pab(«,W n ) = -^-[Vabfpbb — Paa)(u, <J n )] ■ J+(uJ n ,v) 

Paa(U 7 Wn) = |A a («)<5n,0 — ^-[(V ab pba — V b aPab)(«, W n )]| ■ 7+K) U ) (34) 
p hb (v,U n ) = |A b («)5 ni0 — ^-[(VbaPab — V ab pba) (u,^n)] | • J+{uJ n ,v) 

where S n ,o is equal to 1 for n=0 and zero otherwise. The 7 functions are the Lorentzian response of the masing molecules moving at velocity 
v to the radiation field and given as follows: 

J±(LO n ,v) = 1/ jr ± i[w a b — (w a b + Wn) ■ (1 — ^)]| — 1/ |r ± i [ixJab ~ — W„]| 



7±(w n ,w) = 1/ 



T ^ iuj n - (1 



V 



c 



1/[Tt«u„] (35) 



The functions y± (u> n , v) peak at ui n = w a b"u/c, which corresponds to the Doppler shift of the resonance frequency of the masing molecules. 
We note that in arriving at Eqs.[34]we have used the similar approximation as in Section II of Menegozzi & Lamb (1978), i.e. neglecting the 
term involving e -1 "™^ ~ z / c ) w JL p( z , v,cu„). Noting further that d < c,we obtain the final expressions for the Lorentzian response of the 
masing molecules. These expressions are the same as that used in Menegozzi & Lamb (1978). 

In the frequency domaine u„, the interaction term, i.e. the product V • p, can be written as the convolution: 

q— + 00 

V • p(v,LO n ) = U,W n -q) (36) 

q= — 00 

The interaction matrix V between radiation field and the masing molecules in the dipole approximation can be expressed as: 

V = -E ■ d~-^{E{z,t)e-^ b(t - z/c) } -d (37) 

where d is the dipole moment operator for the molecule and in the last step we have made use of the rotating wave approximation. For 
simplicity we assume that the matrix elements of the dipole moment are real, i.e d a b = d a b = d. In the frequency domain the interaction 
term has the form: 

V ab (z,^„) = -~dE(z,u n ) (38) 
Therefore from Eqs.[34]we obtain the expression for the off-diagonal element of the density matrix: 

Pab(",^n) = ^2 ^C^a-q) ■ [Pbb(v,U>q) ~ paa(«,^q)] • 7+ b (^n,u) (39) 
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The equation Eq.l36lbecomes: 

id' 



i d 2 x 

• pib(v,U n ) = — E(lV m +n) ■ E*(0J m -q) ■ [p b b* (v , OI q ) - Paa*(v,(Jq)] ■ J^(u m ,v) (40) 



m . q 

We also obtain similar expression for the complex conjugate term: 



V ab * • p ab (w,tJ n ) = --TT- ^ E'fam-n) ■ Efam-q) ■ [p hh (v,LU q ) - Paa(", Wq)] ■ 7+ i '(w m , 1)) (41) 
m,q 

The population inversion of the maser can be calculated in terms of the elements of the density matrix as follows: 

Ap(v,LU n ) = p m (v,U a ) — p hh ((jj n ,v) 

d 2 f 

= 7^2 { y^-B(t^m+B.)^*(^m-q) [Pbb(^q^) ~ Paa(^q, «)] 7- & (^m, v) 



* m > q 



2ft 2 

+ y^-E*(w m _ n )-E(o; m _q) [pbb(Wq,«) - paa(^q, f)] 7+ 6 (^ m , 1)) + 
m, q 

(A a (v) - A b («))<5n,o| 7+( W n) (42) 

The polarisation vector P(z, i) of the masing medium is the response of the molecules to the presence of the electromagnetic field. Quantum 
mechanically, P(jz, t) is the average value of the dipole moment operator d and can be calculated as the trace of p ■ d: 

+oo 

P(;M) = J dvtr(p-d) (43) 

— oo 

or written explicitly using the rotating wave approximation: 

1 f + °° 

-P(z,t)exp{-iui ab (t ~ z/c)} ~ / dvp ab (z,t)d (44) 



2 

Transforming the above equation into frequency domain and make use of Eq.[39]we obtain: 

id 2 f +a ° 

P(wn) = — dv Yl ^MMsf) ' lt(un,v) (45) 

J-oo q 

If we define the normalized profiles cj>±, which satisfy J_°° dv cj>±(uj n ,v) — 1, as follows: 

<j}±(cO n ,v) = • ~{±(uJn,v) (46) 

7TC 

the radiation transfer equation Eq.l30lbecomes: 

dE{uj n ) _ 2tt 2 d 2 



dz 



/+oo 
dv 2J £(w n _ q ) Ap(u q , U)<?i>+ (W n ,w) (47) 
-oo 



2.3 First order approximation 

The equations derived in the previous section relating the transfer of the radiation field to the density matrix of masing molecules are rather 
complicated. One technique to solve them is the perturbative expansion consisting in solving the equations through a series expansion. To 
the first order, the populations are considered constant over time within a given realization. Therefore the population inversion can be written 
simply as Ap(v) with no dependence on harmonic frequency. The density matrix element connecting the upper level to the ground level 
therefore contain electric field to the first order. The transfer equations are particularly simple. 

d£(w„) _ 2tt 2 d 2 



dz h 



dvE(uj n )Ap(v)cj> + (w n ,v) (48) 



The above equation can be cast into more familiar form involving the intensity of the radiation field per unit frequency v (v = lu/2tt) by 
defining I(u} m )Au)/2ir — c/8ir E(u> m )E* (ui m ) and the real part <^ r (w m , v) of <^±(u> m , v) : 

dl(uj n ) 4-7r 2 d 2 f , ,. , . , , . . 

K dz ' = — — / dvl(u n ) ■ Ap(v) ■ 4> T (LUn,v) (49) 

The equation for the population inversion Eq.[42]in this case can be written in a particularly simple form: 
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FAp(v) = 7 ^ ' Ap( - V "> ' {f + + + AA(u) (50) 

n 

where AA(w) = A a (u) — Xb(v). The coefficient 2d 2 /c?i 2 can be converted to the Einstein coefficient B by the relation 8ir 2 d 2 /h 2 c — 
3B/4n 2 . We can simplify the equation if we define the real part of Lorentz functions ■y±'(u) ul , v) as 7 r (a; m , v): 

TAp(v) = • 7 K) ' M") ' Y(un,v)Aca + AX(v) (51) 

n 

We note that the function 7 r ((jj n , w) has the normalization J_ °° 7 r (oj n , du = 7r. From the radiation transfer equation we can readily 
derive the unsaturated optical depth of the maser as: 

T0 (v) = ±™ L ^1 = • SB • £) L ^ (52) 

h r V47T v } r 

where L is the length of the maser. The appearance of the factor p/c is due to the fact that p(v) has the dimension of particle number per 
unit velocity (cm _3 /cms _1 ) and needs to be converted to particle number per unit frequency (cm _3 /Hz) with the factor v/c in the above 
equation. 

The equations Eqs.[49]&[5T]can be easily seen as identical to the equation derived by Casperson & Yariv (1972) to describe the amplification 
of optical continuum in laser amplifier with both homogeneous and inhomogeneous broadening. 

The standard radiation transfer equation for the average intesity of maser emission can be recovered by taking ensemble-averaging of the 
above equations and postulating that the maser medium is in steady state with fluctuations of Ap(v) independent from that of J(a; m ). 
Of course the last assumption is very drastic because the population inversion is driven by the radiation field. However, Goldreich et al. 
(1973) suggest that in the small signal regime the fluctuations of the population inversion is small and has very narrow bandwidth, or 
the autocorrelation time much longer than that of the radiation field. That fact can justify the assumption of ignoring the fluctuations of 
population inversion with respect to that of the radiation field. Noting further that Y and <f> v are sharply peaked functions and satisfy the 
normalization ^ n 7 r (cj n , v)Alo — tv and f^°° (j> 1 (uj n ,v)dv — 1, we obtain the standard equation describing the ensemble anverage 
intensity < I(u> n > for maser: 

d < I(iOn) > _ r( > Tp(v) 

d{z/L) - <I ^> > i + (3B/ 27 rr) < I(u n ) > (53) 

This equation was also derived by Litvak (1970) using the integral equation approach, which makes use explicitly of the assumption on 
Gaussian statistics of the radiation field. 



2.4 Change of radiation statistics 

To characterize the deviation from Gaussian statistics of the radiation field we define the parameter S following Menegozzi & Lamb (1978): 

r/ \ < /2 ( w m) > - < J(w m ) > 2 

5 Wm = 6 rt wT " < 54 ) 

For complex Gaussian random variable the value of S is equal to unity. The evolution of S in the unsaturated regime can be derived heuristi- 
cally by taking the derivative of 5: 

1 (±<l\^)>-2 ^^ ± <!(.)>) (55) 



dz < I(oj m ) > 2 \dz < /(w m ) > dz 

To estimate the change of statistics of the maser radiation field, we restrict ourselves to the first order approximation as shown in prevous 
section, mainly the equation Eq.|49] In this approximation the population inversion is considered to be constant within a given realization. 
We note further that because the function r (w n , v) is a sharply peaked, the population inversion Ap(v) at any velocity v is mainly affected 
by the Fourier component cu n of the radiation field in resonance with the molecular transition. Thus, we may consider the function (j> r (io n , v) 
as a delta function in velocity and ignore the contribution of molecules having velocities not in resonance with the radiation field. In this 
case, the equations for the ensemble average < I(u) m ) > and < J 2 (tj m ) > follow easily: 

d 47r 2 d 2 

— < > = — — < I(w m ) Ap > (56) 

^ < I 2 K0 > = < / 2 K0 Ap > (57) 

It should be noted here that the fluctuations of Ap are induced by the radiation field. In the unsaturated or small signal regime we can 
approximate the dependence of population inversion on the intensity of the radiation field as: 



Ap ~ — 
1 V 



i _ ^Lnuj m ) 



(58) 

Substituting the above expressions into Eq.[55]we obtain: 
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dS_ ftrMV* AA 3B ( < ^ > + < /> m ) > 2 (59) 



< /(waO > 2 r 27rr y v y < 7(w m ) > 

In the small signal regime, the initial radiation field has Gaussian statistics. As shown in Menegozzi & Lamb (1978), for Gaussian variables 
we can factorize the terms involving power of intensity / such as < / 3 (oj m ) > = 6 < I(uJ m ) > 3 and < I 2 {ujm) > = 2 < I(u) m ) > 2 . 
Therefore the change in statistics of the total radiation field seen by the molecules inside the maser is: 

dS 167r 2 d 2 AA W T , . 

5 = ^-2^r </( ^ )> (60) 

The above equation shows that the parameter 8 will decrease as the radiation is amplified by the maser. That means large fluctuations of 
the field are suppressed in comparision to smaller fluctuations. The result has a simple physical interpretation: large fluctuations can quickly 
deplete the population inversion, resulting in slower growth rate. Smaller fluctuations can grow at a different and faster rate, thus effectively 
reduce the range of intensity fluctuations. 

We can also assess heuristically the effect of the interaction between masing molecules and the radiation field by comparing the intensity 
predicted by the standard radiation transfer equation Eq.[53]and our model in the first order approximation described by Eqs.[56]&[58] The 
amplification of the ensemble averaged intensity of the radiation field in this approximation is: 

d r , , 47r 2 d 2 AA / T , . W r2 . . 

The appearance of factor 2 in the second term indicates that the fluctuations of radiation field can reduce the population inversion faster then 
that predicted by the standard equation Eq.[53] Consequently, our simulations will give a slightly lower ensemble average intensity under the 
same conditions for the maser in comparision to the standard maser theory. We note that previous work by Field & Richardson (1984) and 
Field & Gray (1988), taking into account in an approximate manner the partial coherence of the maser radiation, also reached qualitatively 
similar conclusion, namely lower maser intensity appears in modest saturation regime in comparison to the standard maser theory. 



3 SIMULATION OF RADIATION TRANSPORT 

In our simulation we will choose the parameters appropriate to astronomical maser. For simplicity, we choose the loss rate T = 1 s _1 and 
the normalized pump rates are assumed have a velocity dispersion a with the difference in pump rates AX(v) = exp(— v 2 /a 2 ). In our 
simulation we use 200 modes with a frequency resolution of Alo = 0.75 F or 0.75 s _1 around the maser line and we assume that this 
frequency range covers the range — a to +a in the velocity domain. The actual velocity resolution Av will depend on the frequency of the 
maser line, as Av = c/u) a bAw. For the maser line such as the 1612 MHz OH maser, Av is approximately 2 cm s _1 . The corresponding 
value for the velocity dispersion is 200 cm s _1 Although the velocity dispersion of the maser line in our simulations is much smaller than 
in real astronomical masers, the number of frequency modes and the bandwidth are large enough, i.e. the bandwidth of 150 s -1 is much 
greater than the loss rate of 1 s _1 , to capture the main features of the broadband radiation field produced by the astronomical masers. Since 
we deals with only partially saturated maser, we consider only 15 harmonic components of the population inversion p(u}„)- As shown later, 
the number of harmonic components is enough to capture the pulsation of the molecular population inversion. We generate the background 
continuum radiation in a simimlar way as Menegozzi & Lamb (1978) using random generator ran2 from Press et al. (1992). The phase 
of electric field components is random and uniformly distributed over the interval to 2n. For the sake of simplicity, we assume that the 
amplitude of different modes has a Gaussian distribution with zero mean and variance < I(uj n )Au) > = 1 on scale of 2hu 3 /c 2 . On this scale 
the intensity has the unit of the photon occupation number. The intensity of the maser shown in all figures is also of the form I(ui n )Au, 
where I(oj n ) is the intensity per unit frequency v as defined in previous section. We use fourth-order Runge-Kutta method with a fixed step 
h = 0.01 to integrate the Eq.[47] Each realization is evolved through the maser using the rate equations Eqs.[42]together with the radiation 
transfer equation Eq.|47]and we record the emergent radiation for later analysis. 

We carry out our simulations in the partially saturated regime, which is likely relevant to most astronomical masers. The choice is also 
necessary because we consider only a limited number of harmonic components of the molecular population inversion. In the saturated 
regime the fluctuation will be stronger and thus require the consideration of a larger number of harmonic components. The amplification of 
the maser is specified by the unsaturated optical depth ro(v — 0) at the line center. A value of tq(v — 0) = 22 is used throughout in our 
simulations. The spontaneous transition rate between the upper and lower maser levels is taken to be 10~ 9 s _1 . The Einstein coefficient B is 
related to the spontaneous transition rate by the well-known relation A = 2hv 3 /c 2 B. 



4 DISCUSSION 

The new ingredient of our model is the explicit treatment of the radom radiation field and the inclusion of the molecular population pulsation 
and the coupling between different modes of the radiation field. We can examine each of these effects on the propagation of electromagnetic 
waves in the maser. 
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If we neglect the population pulsation, as can be seen from Eq. [47] only the intensity of the waves enters into the radiation transfer 
equation. Thus different realizations with initial constant amplitude but random phases are indistinguishable and we expect them to evolve 
similarly. However, when population pulsation is taken into account, different modes can interact and the random phases can induce fluc- 
tuations in the amplitude of the waves. In Fig.Q we show the evolution of a particular realization of the field whose modes have constant 
amplitude but random phases. Initially the mode-coupling is negligible and the intensity across the maser line half way into the maser 
medium, which corresponds to the optical depth of r/2, varies smoothly according to the change of the optical depth with frequency. The 
usual line-narrowing effect can be clearly seen because the initial FWHM of the maser line due to inhomogeneous (Doppler) broadening 
is 2<rln2~ 140 frequency channels and the FWHM of the laser line at the mid-point of the laser (see the middle frame of Fig. QJ is only 
around 50 frequency channels. As the intensity grows different modes of the radiation field inside the homogeneous Lorentzian profile can 
couple and drive the molecular population pulsations. The emergent maser radiation fluctuates widely around the line center where this effect 
is strongest (Fig. |T|l. The maser in this case is only partially saturated with the saturation parameter R/F ~ 0.8, where R is the simulated 
emission rate. In Fig. [2] we show the amplitude of the Fourier components for the population pulsations at the line center. The dominant 
component is at zero frequency or the average of population inversion. The amplitude of other Fourier components remains small and de- 
creases at higher pulsation frequency. Nevertheiiess the population pulsation and mode-coupling are effective in inducing the fluctuations of 
the amplified maser radiation. We can also see from Fig.Q]that the phases of different modes remain random as the waves propagate through 
the maser and experience amplification. No mode-locking or phase correlation is seen in our simulations. We also checked the simulation 
results shown in Menegozzi & Lamb (1978) and could not find any evidence for phase correlation between different modes of the radiation 
field. Thus, our results do not corroborate the suggestion of Gray & Bewley (2003) that there is a strong phase correlation between frequency 
modes or mode-locking even for partially saturated maser. The likely reason for the difference is that the formulation of Gray & Bewley 
(2003), which follows closely that of Yariv (1989) for the amplification of monochromatic (single mode) laser emission, instead of the broad- 
band random radiation field. The steady state constraint imposed on Eq. 28 of their paper leads to the presence in subsequent equations of 
time-dependent terms, which are supposed to be eliminated by the same procedure. Of course, the problem can be avoided if their equation 
Eq. 28 is expanded into Fourier series as we have done here. 

In Fig. [3] we show the input radiation field with both random amplitude and phase of different modes are random. We also show in Fig. [3] 
the radiation field at the mid-point of the maser and the emergent radiation field. As can be seen in the figure, only modes around the line 
center are amplified strongly. Because of the nonlinear amplification, the relative intensity fluctuations between different mode changes as 
the radiation is amplified in the maser medium. That is the case for the relatively strong mode seen close to the line center in the initial input 
field (see Fig.[3]>. When the intensity of the radiation field is still small, i.e. the maser is unsaturated, all the modes are amplified by the same 
factor as evident in the middle frame of Fig.[3] However, as the maser approaches saturation, the adjacent modes, which are initially weaker, 
grow faster. As a result, the relative intensity fluctuation is reduced in the output radiation field. We also note that the phases between modes 
are still completely random as seen in Fig. [5] 

We have followed the amplification of an ensemble of 1 8000 realizations of the incident radiation field through the maser. By taking the 
ensemble average of the output field we could build up the observable intensity profile of the maser line and the statistics of the radiation field. 
The parameter 5, which measures the deviation from Gaussian statistics, is shown in Fig. [4] For the first time we are able to determine the 
statistics of the maser radiation directly from numerical simulations. Even in our case of partial saturation with R/r ~ 0.8 at the line center, 
the statistical distribution of the total radiation field deviates quite significantly from Gaussian. The decrease of S at the line center suggests 
that intensity fluctuations diminish with saturation and the statistical distribution of the radiation field becomes more centrally peaked. As 
discussed previously, the non-Gaussian statistics is the direct consequence of the non-linear interaction between masing molecules and the 
radiation field. We could expect that the deviation will be even more pronounced as the maser becomes fully saturated. 
We also compare the ensemble averaged intensity of the total electric field taken over 18000 realizations of the field (see Fig. |4]> with the 
prediction of the standard radiation transfer equations Eq.|53] The two calculations match closely with each other in the line wings. At the 
center of the maser line our simulation produces a slightly lower ensemble average intensity than that given by the standard theory, consistent 
with our heuristic analysis in previous section. Therefore in the unsaturated and partially saturated regime the standard theory of astronomical 
maser can produce results in good agreement with that from our more realistic treatment. However, in saturated masers where population 
pulsation and mode-coupling are expected to be much stronger than considered here and the radiation field deviates even more from Gaussian 
statistics, the usual approach used by Litvak (1970) and Goldreich et al. (1973) are no longer applicable. As a result, previous calculations 
using the standard radiation transfer equations and carried out in the full saturation regime with R/T 2> 10 (Western & Watson 1984, 
Deguchi & Watson 1990, Nedoluha & Watson 1990) might not give accurate description of the maser radiation and should be considered 
with cautions. We also note that fully saturated maser is of purely theoretical interest because real astronomical masers are unlikely to be 
found in this regime due to various constraints including the pumping of masers. Although it is difficult to determine the degree of saturation 
for astronomical masers, the modelling results of Watson et al. (2002) suggest that recent observations of water maser features with Gaussian 
spectral line profile (Sarma et al. 2001) can only be explained if the maser is unsaturated with R/r < 1/3. 

The deviation from Gaussian statistics seen in our simulations can be measured in practice if we can observe an isolated maser source with 
frequency resolution comparable to the homogeneous line width ~ F. When the frequency resolution is much larger than the homogenous 
bandwidth, the measured electric field is the sum over many independently random harmonic components and the resulting statistics is ex- 
pected to become Gaussian. Large bandwidth is used in previous attempts (Evans et al. 1972, Moran 1981) to measure the statistics of the 
maser radiation and that might contribute to the negative results, although other complications such as the presence of several independent 
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maser soures inside the radio telescope beam might preclude any such detection. 

Finally we should emphasize that our formulation of the maser theory, which is based on the work of Menegozzi & Lamb (178), as 
well as the standard theory of maser of Litvak (1970) and Goldreich et al. (1973) are all based on the one-dimensional idealization of the 
maser medium. Thus the maser radiation field consists of plane waves propagating along the same direction. In real astronomical masers, the 
masing medium is a three-dimensional object. The radiation field at any point in the medium is the superposition of waves propagating in 
many different directions. The formaluation of the maser theory in three dimensions, although similar in principle, will be more complicated. 
The abovementioned results and conclusions from our simulations will need to be reconfirmed in the future with fully three dimensional 
calculations. In addition, the spontaneous emission by the masing molecules has been omitted in our model. This simplification is quite 
common in the theoretical studies of astronomical masers. An explicit treatment of the spontaneous emission will require the reformulation 
of the maser theory using quantum electrodynamics, probably along the line presented in the classic work of Sargent et al. (1974) and clearly 
is out of the scope of our paper. However, we note that, in a number of cases the astronomical masers are known observationally to amplified 
the background continuum source (such as the central red giant star or the AGN). Thus it's reasonable to expect that the omission of the 
spontaneous emission will not change qualitatively the results of our model. 



5 CONCLUSION 

The explicit incorporation of the broadband random radiation field and the molecular population pulsation into our treatment of astronomical 
maser has proved very important to investigate the properties of maser emission. We have shown clearly the effect of mode coupling in 
driving the molecular population pulsation. The amplification of background radiation by maser is accompanied by the change of statistics 
of the radiation field in which the large fluctuations are suppressed relative to smaller fluctuations. Our simulation results suggest that 
the standard radiation transfer equation provides a good description of the maser properties such as intensity and line-narrowing in the 
unsaturated and partially saturated regime. Howerver, the application of the same equation to study masers in the strong saturation regime 
should be considered with caution. Further studies on the effect of mode coupling and population pulsation in saturated masers are needed 
and we hope to address them in a future publication. 
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Figure 1. Evolution of a realization of the radiation field with different modes having the same amplitude but random phases. Intensity is on the scale of 
2hv 3 /c 2 . Left frame shows the initial intensity and phases of the radiation field. Middle frame shows the intensity and phases of the field half-way into the 
maser medium, corresponding to To (v = 0) = 1 1. Right frame shows the intensity and phases of radiation field at the output of the maser. 
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Figure 2. Amplitude of Fourier components of the population pulsation with velocity v = at the output of the maser medium (z = L), corresponding to the 
simulation presented in Fig.[T] 
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Figure 3. Evolution of a realization of the radiation field with different modes having the random amplitude and random phases. Intensity is on the scale of 
2hv /c . Left frame shows the initial intensity and phases of the radiation field. Middle frame shows the intensity and phases of the field half-way into the 
maser medium, corresponding to tq (v = 0) = 1 1. Right frame shows the intensity and phases of the emerging radiation field at the output of the maser. 
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Figure 4. Upper frame: ensemble average power spectrum taken over 18000 independent samples of the emerging radiation field. The maser emission 
calculated using the standard radiation transfer equation is shown in dotted line. Intensity is on the scale of 2hv 3 /c 2 . Lower frame: ensemble average of 
parameter 5, i.e. the departure from Gaussian statistics, taken over the same number of independent samples of the emerging radiation field. 



